﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides the numerical integration of infinity integrals by using the Gauss-Kronod method.
    /// </summary>
    [Serializable]
    public class GaussKronrodInfinityIntegrator
    {
        /// <summary>
        /// The integral for the numerical integration.
        /// </summary>
        private InfinityDoubleIntegral integral;

        /// <summary>
        /// Saves the absolute error of the numerical integration.
        /// </summary>
        private double absoluteError;

        /// <summary>
        /// Saves the absolute accuracy of the numerical integration.
        /// </summary>
        private double absoluteAccuracy;

        /// <summary>
        /// Saves the relative accuracy of the numerical integration.
        /// </summary>
        private double relativeAccuracy;

        /// <summary>
        /// Describes if the numerical integration has a precision error.
        /// </summary>
        private bool precisionError;

        /// <summary>
        /// Saves the needed function evaluations of the numerical integration.
        /// </summary>
        private double neededFunctionEvaluations;

        /// <summary>
        /// Initializes a new instance of the <see cref="GaussKronrodInfinityIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The integral for the numerical integration.</param>
        public GaussKronrodInfinityIntegrator(InfinityDoubleIntegral integral)
        {
            this.integral = integral;
        }

        /// <summary>
        /// Gets or sets the integral for the numerical integration.
        /// </summary>
        /// <value>The integral for the numerical integration.</value>
        public InfinityDoubleIntegral Integral
        {
            get { return integral; }
            set { integral = value; }
        }

        /// <summary>
        /// Gets the absolute error of the numerical integration.
        /// </summary>
        /// <value>The absolute error of the numerical integration.</value>
        public double AbsoluteError
        {
            get { return absoluteError; }
        }

        /// <summary>
        /// Gets the absolute accuracy of the numerical integration.
        /// </summary>
        /// <value>The absolute accuracy of the numerical integration.</value>
        public double AbsoluteAccuracy
        {
            get { return absoluteAccuracy; }
        }

        /// <summary>
        /// Gets the relative accuracy of the numerical integration.
        /// </summary>
        /// <value>The relative accuracy of the numerical integration.</value>
        public double RelativeAccuracy
        {
            get { return relativeAccuracy; }
        }

        /// <summary>
        /// Gets the needed function evaluations of the numerical integration.
        /// </summary>
        /// <value>The needed function evaluations of the numerical integration.</value>
        public double NeededFunctionEvaluations
        {
            get { return neededFunctionEvaluations; }
        }

        /// <summary>
        /// Describes if the numerical integration has a precision error. If true the
        /// calculated result may be not correct.
        /// </summary>
        /// <value><c>true</c> if [precision error]; otherwise, <c>false</c>.</value>
        public bool PrecisionError
        {
            get { return precisionError; }
        }

        /// <summary>
        /// Integrates an infinity integral by using the Gauss-Kronrod method.
        /// </summary>
        /// <param name="epsabs">The absolute requested accuracy of the result.</param>
        /// <param name="epsrel">The relative requested accuracy of the result.</param>
        /// <returns>The area of the specified integral.</returns>
        public double Integrate(double epsabs, double epsrel)
        {
            int id = 0;
            int k;
            int inf = 0;
            int ier = 0;
            int neval = 0;
            int jupbnd = 0;
            const int interations = 500;
            int limit = interations - 1;
            double abserr = 0;
            double bound = this.integral.Bound;
            double small = 0;
            double abseps = 0;
            double area1 = 0;
            double area12 = 0;
            double area2 = 0;
            double a1 = 0;
            double a2 = 0;
            double b1 = 0;
            double b2 = 0;
            double correc = 0;
            double defabs = 0;
            double defab1 = 0;
            double defab2 = 0;
            double erlarg = 0;
            double erlast = 0;
            double error1 = 0;
            double error2 = 0;
            double erro12 = 0;
            double errsum = 0;
            double ertest = 0;
            double resabs = 0;
            double reseps = 0;
            double[] res3la = new double[3];
            double[] alist = new double[interations];
            double[] blist = new double[interations];
            double[] elist = new double[interations];
            double[] rlist = new double[interations];
            double[] rlist2 = new double[52];
            int[] iord = new int[interations];

            switch (this.integral.Range)
            {
                case IntegralRange.MinusInfinityToBound:
                    inf = -1;
                    break;

                case IntegralRange.BoundToPlusInfinity:
                    inf = 1;
                    break;

                case IntegralRange.MinusInfinityToPlusInfinity:
                    inf = 2;
                    break;
            }

            ier = 0;
            neval = 0;
            int last = 0;
            double result = 0.0;
            abserr = 0.0;
            alist[0] = 0.0;
            blist[0] = 1.0;
            rlist[0] = 0.0;
            elist[0] = 0.0;
            iord[0] = 0;
            if ((epsabs < 0.0) && (epsrel < 0.0))
            {
                ier = 6;
            }
            if ((inf != 1) && (inf != -1) && (inf != 2))
            {
                ier = 6;
            }
            if (ier == 6)
            {
                return result;
            }

            double boun = bound;
            if (inf == 2)
            {
                boun = 0.0;
            }

            result = GaussKronrodRules.G_K15I(this.integral.Function, boun, inf, 0.0, 1.0, ref abserr, ref defabs,
                                              ref resabs);
            last = 0;
            rlist[0] = result;
            elist[0] = abserr;
            iord[0] = 0;
            double dres = Math.Abs(result);
            double errbnd = (((epsabs) > (epsrel * dres)) ? (epsabs) : (epsrel * dres));
            if ((abserr <= 100.0 * Double.Epsilon * defabs) && (abserr > errbnd))
            {
                ier = 2;
            }
            if (limit == 0)
            {
                ier = 1;
            }
            if ((ier != 0) || (abserr <= errbnd) && (abserr != resabs) || (abserr == 0.0))
            {
                goto _130;
            }

            rlist2[0] = result;
            double errmax = abserr;
            int maxerr = 0;
            double area = result;
            errsum = abserr;
            abserr = Double.MaxValue;
            int nrmax = 0;
            int nres = 0;
            int ktmin = 0;
            int numrl2 = 1;
            int extrap = 0;
            int noext = 0;
            int ierro = 0;
            int iroff1 = 0;
            int iroff2 = 0;
            int iroff3 = 0;
            int ksgn = -1;
            if (dres > (1.0 - 50.0 * Double.Epsilon) * defabs)
            {
                ksgn = 1;
            }

            for (last = 1; last <= limit; last++)
            {
                a1 = alist[maxerr];
                b1 = 0.5 * (alist[maxerr] + blist[maxerr]);
                a2 = b1;
                b2 = blist[maxerr];
                erlast = errmax;
                area1 = GaussKronrodRules.G_K15I(this.integral.Function, boun, inf, a1, b1, ref error1, ref resabs,
                                                 ref defab1);
                area2 = GaussKronrodRules.G_K15I(this.integral.Function, boun, inf, a2, b2, ref error2, ref resabs,
                                                 ref defab2);
                area12 = area1 + area2;
                erro12 = error1 + error2;
                errsum = errsum + erro12 - errmax;
                area = area + area12 - rlist[maxerr];
                if ((defab1 == error1) || (defab2 == error2))
                {
                    goto _15;
                }
                if ((Math.Abs(rlist[maxerr] - area12) > 1.0e-5 * Math.Abs(area12)) || (erro12 < .99 * errmax))
                {
                    goto _10;
                }
                if (extrap != 0)
                {
                    iroff2++;
                }
                else
                {
                    iroff1++;
                }
                _10:
                if ((last > 9) && (erro12 > errmax)) // last > 10 
                {
                    iroff3++;
                }
                _15:
                rlist[maxerr] = area1;
                rlist[last] = area2;
                errbnd = (((epsabs) > (epsrel * Math.Abs(area))) ? (epsabs) : (epsrel * Math.Abs(area)));

                if (((iroff1 + iroff2) >= 10) || (iroff3 >= 20))
                {
                    ier = 2;
                }
                if (iroff2 > 5)
                {
                    ier = 3;
                }

                if (last == limit)
                {
                    ier = 1;
                }

                if ((((Math.Abs(a1)) > (Math.Abs(b2))) ? (Math.Abs(a1)) : (Math.Abs(b2))) <=
                    (1.0 + 1000.0 * Double.Epsilon) * (Math.Abs(a2) + 1000.0 * Double.MinValue))
                {
                    ier = 4;
                }

                if (error2 <= error1)
                {
                    alist[last] = a2;
                    blist[maxerr] = b1;
                    blist[last] = b2;
                    elist[maxerr] = error1;
                    elist[last] = error2;
                }
                else
                {
                    alist[maxerr] = a2;
                    alist[last] = a1;
                    blist[last] = b1;
                    rlist[maxerr] = area2;
                    rlist[last] = area1;
                    elist[maxerr] = error2;
                    elist[last] = error1;
                }

                DqSort(limit, last, ref maxerr, ref errmax, elist, iord, ref nrmax);
                if (errsum <= errbnd)
                {
                    goto _115;
                }
                if (ier != 0)
                {
                    goto _100;
                }
                if (last == 1)
                {
                    goto _80;
                }
                if (noext != 0)
                {
                    continue;
                }
                erlarg -= erlast;
                if (Math.Abs(b1 - a1) > small)
                {
                    erlarg += erro12;
                }
                if (extrap != 0)
                {
                    goto _40;
                }

                if ((Math.Abs(blist[maxerr] - alist[maxerr])) > small)
                {
                    goto _90;
                }
                extrap = 1;
                nrmax = 1;
                _40:
                if ((ierro == 3) || (erlarg <= ertest))
                {
                    goto _60;
                }

                id = nrmax;
                jupbnd = last;
                if (last > (2 + limit / 2))
                {
                    jupbnd = limit + 3 - last;
                }
                for (k = id; k <= jupbnd; k++)
                {
                    maxerr = iord[nrmax];
                    errmax = elist[maxerr];
                    if (Math.Abs(blist[maxerr] - alist[maxerr]) > small)
                    {
                        goto _90;
                    }
                    nrmax++;
                }

                _60:
                numrl2++;
                rlist2[numrl2] = area;
                reseps = DqExt(ref numrl2, rlist2, ref abseps, res3la, ref nres);
                ktmin++;
                if ((ktmin > 5) && (abserr < 1.0e-3 * errsum))
                {
                    ier = 5;
                }
                if (abseps >= abserr)
                {
                    goto _70;
                }
                ktmin = 0;
                abserr = abseps;
                result = reseps;
                correc = erlarg;
                ertest = (((epsabs) > (epsrel * Math.Abs(reseps))) ? (epsabs) : (epsrel * Math.Abs(reseps)));
                if (abserr <= ertest)
                {
                    goto _100;
                }

                _70:
                if (numrl2 == 0)
                {
                    noext = 1;
                }
                if (ier == 5)
                {
                    goto _100;
                }
                maxerr = iord[0];
                errmax = elist[maxerr];
                nrmax = 0;
                extrap = 0;
                small = small * 0.5;
                erlarg = errsum;
                continue;
                _80:
                small = .375;
                erlarg = errsum;
                ertest = errbnd;
                rlist2[0] = area;
                _90:
                ;
            }
            _100:
            if (abserr == Double.MaxValue)
            {
                goto _115;
            }
            if ((ier + ierro) == 0)
            {
                goto _110;
            }
            if (ierro == 3)
            {
                abserr += correc;
            }
            if (ier == 0)
            {
                ier = 3;
            }
            if ((result != 0.0) && (area != 0.0))
            {
                goto _105;
            }
            if (abserr > errsum)
            {
                goto _115;
            }
            if (area == 0.0)
            {
                goto _130;
            }
            goto _110;
            _105:
            if (abserr / Math.Abs(result) > errsum / Math.Abs(area))
            {
                goto _115;
            }

            _110:
            if ((ksgn == -1) &&
                ((((Math.Abs(result)) > (Math.Abs(area))) ? (Math.Abs(result)) : (Math.Abs(area))) <= defabs * .01))
            {
                goto _130;
            }
            if ((0.01 > result / area) || (result / area > 100.0) || (errsum > Math.Abs(area)))
            {
                ier = 6;
            }
            goto _130;

            _115:
            result = 0.0;
            for (k = 0; k <= last; k++)
            {
                result += rlist[k];
            }
            abserr = errsum;
            _130:
            neval = 30 * last + 15;
            if (inf == 2)
            {
                neval *= 2;
            }
            if (ier > 2)
            {
                (ier)--;
            }

            this.absoluteAccuracy = epsabs;
            this.relativeAccuracy = epsrel;
            this.absoluteError = abserr;
            this.neededFunctionEvaluations = neval;
            this.precisionError = !(ier == 0);

            return result;
        }

        private static void DqSort(int limit, int last, ref int maxerr, ref double ermax, double[] elist, int[] iord,
                                   ref int nrmax)
        {
            int i;
            int j;
            int isucc;

            if (last > 1)
            {
                goto _10;
            }

            iord[0] = 0;
            iord[1] = 1;
            goto _90;

            _10:
            double errmax = elist[maxerr];

            if (nrmax == 0)
            {
                goto _30;
            }

            int ido = (nrmax) - 1;

            for (i = 0; i <= ido; i++)
            {
                isucc = iord[nrmax - 1];

                if (errmax <= elist[isucc])
                {
                    goto _30;
                }

                iord[nrmax] = isucc;
                (nrmax)--;
            }

            _30:
            int jupbn = last;

            if (last > (limit / 2 + 2))
            {
                jupbn = limit + 3 - last;
            }

            double errmin = elist[last];
            int jbnd = jupbn - 1;
            int ibeg = nrmax + 1;

            if (ibeg > jbnd)
            {
                goto _50;
            }

            for (i = ibeg; i <= jbnd; i++)
            {
                isucc = iord[i];

                if (errmax >= elist[isucc])
                {
                    goto _60;
                }

                iord[i - 1] = isucc;
            }

            _50:
            iord[jbnd] = maxerr;
            iord[jupbn] = last;
            maxerr = iord[nrmax];
            ermax = elist[maxerr];

            return;

            _60:
            iord[i - 1] = maxerr;

            int k = jbnd;

            for (j = i; j <= jbnd; j++)
            {
                isucc = iord[k];
                if (errmin < elist[isucc])
                {
                    iord[k + 1] = last;
                    maxerr = iord[nrmax];
                    ermax = elist[maxerr];

                    return;
                }

                iord[k + 1] = isucc;
                k--;
            }

            iord[i] = last;

            _90:
            maxerr = iord[nrmax];
            ermax = elist[maxerr];

            return;
        }

        private static double DqExt(ref int n, double[] epstab, ref double abserr, double[] res3la, ref int nres)
        {
            double delta1;
            double delta2;
            double delta3;
            double epsinf;
            double error;
            double err1;
            double err2;
            double err3;
            double e0;
            double e1;
            double e1abs;
            double e2;
            double e3;
            double res;
            double ss;
            double tol1;
            double tol2;
            double tol3;
            int i;
            int ib2;
            int ie;
            int indx;
            int k1;
            int k2;
            int k3;

            (nres)++;
            int NN = n;
            NN++;
            abserr = Double.MaxValue;
            double result = epstab[n];
            if (NN < 3)
            {
                goto _100;
            }
            int limexp = 50;
            epstab[n + 2] = epstab[n];
            int newelm = (n) / 2;
            epstab[n] = Double.MaxValue;
            int num = NN;
            k1 = NN;
            for (i = 1; i <= newelm; i++)
            {
                k2 = k1 - 1;
                k3 = k1 - 2;
                res = epstab[k1 + 1];
                e0 = epstab[k3 - 1];
                e1 = epstab[k2 - 1];
                e2 = res;
                e1abs = Math.Abs(e1);
                delta2 = e2 - e1;
                err2 = Math.Abs(delta2);
                tol2 = (((Math.Abs(e2)) > (e1abs)) ? (Math.Abs(e2)) : (e1abs)) * Double.Epsilon;
                delta3 = e1 - e0;
                err3 = Math.Abs(delta3);
                tol3 = (((e1abs) > (Math.Abs(e0))) ? (e1abs) : (Math.Abs(e0))) * Double.Epsilon;
                if ((err2 > tol2) || (err3 > tol3))
                {
                    goto _10;
                }
                result = res;
                abserr = err2 + err3;
                goto _100;
                _10:
                e3 = epstab[k1 - 1];
                epstab[k1 - 1] = e1;
                delta1 = e1 - e3;
                err1 = Math.Abs(delta1);
                tol1 = (((e1abs) > (Math.Abs(e3))) ? (e1abs) : (Math.Abs(e3))) * Double.Epsilon;
                if ((err1 <= tol1) || (err2 <= tol2) || (err3 <= tol3))
                {
                    goto _20;
                }
                ss = 1.0 / delta1 + 1.0 / delta2 - 1.0 / delta3;
                epsinf = Math.Abs(ss * e1);
                if (epsinf > 1.0e-4)
                {
                    goto _30;
                }
                _20:
                NN = i + i - 1;
                goto _50;
                _30:
                res = e1 + 1.0 / ss;
                epstab[k1 - 1] = res;
                k1 -= 2;
                error = err2 + Math.Abs(res - e2) + err3;
                if (error > (abserr))
                {
                    goto _40;
                }
                abserr = error;
                result = res;
                _40:
                ;
            }
            _50:
            if (NN == limexp)
            {
                NN = 2 * (limexp / 2) - 1;
            }
            int ib = 1;
            if (((num / 2) * 2) == num)
            {
                ib = 2;
            }
            ie = newelm + 1;
            for (i = 1; i <= ie; i++)
            {
                ib2 = ib + 2;
                epstab[ib - 1] = epstab[ib2 - 1];
                ib = ib2;
            }
            if (num == NN)
            {
                goto _80;
            }
            indx = num - NN + 1;
            for (i = 1; i <= NN; i++)
            {
                epstab[i - 1] = epstab[indx - 1];
                indx++;
            }
            _80:
            if (nres > 3)
            {
                goto _90;
            }
            res3la[(nres) - 1] = result;
            abserr = Double.MaxValue;
            goto _100;
            _90:
            abserr = Math.Abs(result - res3la[2]) + Math.Abs(result - res3la[1]) + Math.Abs(result - res3la[0]);
            res3la[0] = res3la[1];
            res3la[1] = res3la[2];
            res3la[2] = result;
            _100:
            abserr = (((abserr) > (5.0 * Double.Epsilon * Math.Abs(result)))
                          ? (abserr)
                          : (5.0 * Double.Epsilon * Math.Abs(result)));
            n = NN - 1;
            return result;
        }
    }
}